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; ABSTRACT 

" ■ Stellar feedback in galactic bulges plays an essential role in shaping the evolution of 
(N ■ 

galaxies. To quantify this role and facilitate comparisons with X-ray observations, we 

, conduct three-dimensional (3D) hydro dynamical simulations with the adaptive mesh 

O: refinement code, FLASH, to investigate the physical properties of hot gas inside a 

'q^I galactic bulge, similar to that of our Galaxy or M31. We assume that the dynamical 

Q . and thermal properties of the hot gas are dominated by mechanical energy input from 

• supernovae (SNe), primarily Type la, and mass injection from evolved stars as well as 

liHjl iron enrichment from SNe. We study the bulge- wide outflow as well as the SN heating 

7—{ • on scales down to 4 pc. An embedding scheme that is devised to plant individual 

\jQ ' SN remnant (SNR) seeds, allows to examine, for the first time, the effect of sporadic 

00 '. 

, SNe on the density, temperature, and iron ejecta distribution of the hot gas as well 

' as the resultant X-ray morphology and spectrum. We find that the SNe produce a 

(N ■ 

' bulge wind with highly filamentary density structures and patchy ejecta. Compared 

o\'. 

. with a one-dimensional (ID) spherical wind model, the non- uniformity of simulated 

■ gas density, temperature, and metallicity substantially alters the spectral shape and 

^ I increases the diffuse X-ray luminosity. The differential emission measure as a function 

of temperature of the simulated gas exhibits a log-normal distribution, with a peak value 



much lower than that of the corresponding ID model. The X-ray luminosity depends 
sensitively on the mass loss rate from evolved stars. The bulk of the X-ray emission 
comes from the relatively low temperature and low abundance gas shells associated with 
SN blastwaves. SN ejecta are not well mixed with the ambient medium, at least in the 
bulge region. These results, at least partly, account for the apparent lack of evidence for 
iron enrichment in the soft X-ray-emitting gas in galactic bulges and intermediate-mass 
elliptical galaxies. The bulge wind helps to explain the "missing" stellar feedback in 
such galaxies. But the resultant diffuse emission is more than one order of magnitude 
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less than that observed in the Galactic and M31 bulges, indicating that gas in these 
bulges is in a subsonic outflow state probably due to additional mass loading to the 
hot gas and/or due to energy input rate that is substantially lower than the current 
estimate. 

Subject headings: Galaxy: bulge — hydrodynamics — ISM: abundance — ISM: struc- 
ture — supernova remnants: kinematics and dynamics 



Introduction 



Bulges of early-type spirals and elliptical galaxies comprise primarily old low-mass stars, which 
account for more than half of the total stellar mass in the local Universe (Fukugita, Hogan, & Peebles 
1998). These stars collectively generate a long-lasting feedback via mass injection from evolved stars 
(mainly red giant branch stars and planetary nebulae) in the form of stellar winds and energy input 
from Type la SNe. Because of the SN heating, the ISM shoul d be mostly in X-ray-emi t ting plasma 
inside galactic bulges, where little cool gas is present (e.g., iMathews &: Bakerl Il97ll : ISage et al. 



20071 ) . Observations have shown that the X-ray- i nferred gas mass and en ergy are far less than 



those empirical predictions (e.g.. 



Sato &: Tawara 



199S 



. !)avid 



et al 



20061 ). In other words, the 



20071 ). This "missing" stellar feedback 



bulk of stellar feedback expected is not observed (jWang 
problem becomes particularly acute in so-called low Lx/Lb (i.e., the ratio of X-ray luminosity to 
blue band luminosity) bulge-dominated galaxies (typically Sa spirals, SO, and low mass ellipticals). 
After removing the contribution from point sources in those relatively deep Chandra observations, 
the rer naining "diffuse" X-ray component genera lly shows a soft spectrum, indicating a thermal 



origin (jlrwin et al. 



20021 : 



O' Sullivan 



expected Type la SNe energy input ()Li &: Wang 



et al. 



2003 



and its luminosity is only a few percent of the 



2OO7I : 



Li et al 



20071 ). The inferred total mass of 



the X-r ay-emitting gas fall s far short of what is deduced from the stellar mass loss over the galaxy's 



lifetime (David et al 



20061). 



The presence of a bulge-wide outflow may solve the "missing" stellar feedback problem (e.g. 
Tang et al. II2008I ). The ID solution of an SN-heated bulge outflow, however, has problems. The 
physical state of the gas outflow mainly depends on the mass and energy input rates as well as 
the gravitational field. Within those low Lx/Lb galactic bulges of typical mass and energy input 



1991 



.1 



rates , a hot bulge winqj should be present theoretically (e.g., 



Mathews &: Baker 



1971 



Ciotti et al. 



). However, observations of the diffuse hot ISM are apparently at odds with the theoretical wind 
models in nearly all aspects. Firstly, the predicted X-r ay luminosity in a bulge wind scenario is a 
few orders of magnitude smaller than the observed (e.g., iTang et al. II2008I ). Secondly, the expected 



wind temperature is ab out ~1 keV or higher, while the obse rvation-inferred gas temperatures are 
substantially lower (e.g., 



David et al. 


2006 


Li & Wane: 


2007) 



brightness profile of the wind should be steeper than that of starlight (jCiotti et al 



1991 



), but the 

observed profiles of diffus e emission distribut ions are fairly extende d in most low Lx/Lb bulges 
or ellipticals (e.g. fig. 7 in lSarazin et al.ll200ll : fig. 5 in lLi et al.l 120071 ). Furthermore, the predicted 



mean iron abundance of the diffuse hot gas is 3 -7 times solar because of the Type la SN enrichment 



(e.g.; 



Ciotti et al. 



1991 



Sato &: Tawara 



19991 ). whereas the observed spectra usually indicate near- 



or sub-solar iron abundance. 



Some, if not all, of the listed discrepancies may arise from various oversimplifications of the 



exist i ng ID models (e.g., 



Mathews &: Baker 



1987 



1971 



White &: Chevalier 



Ciotti et al. 



1991 



19831 : 



Lowenstein &: Mathews 



). In such models the mechanical energy input of SNe is always treated 
as pure thermal energy smoothly injected into the ISM. In reality, however, SNe, sporadic in both 
time and space, should naturally produce inhomogeneity in the ISM. The density and temperature 
inhomogeneity may significantly affect the X-ray spectrum and luminosity, which are proportional 
to the density square. Explosive energy injection in a hot tenu ous medium can be tr ansported away 
in form of sound wave, so the SN heating is not local (e.g., iTang fc Wang l2005l ) . Furthermore, 
whether or not the SN ejecta of each individual SNR can be well mixed with the surrounding 
material is crucial to address the apparent low abundance puzzle of the hot gas. These effects need 
to be quantified in order to correctly interpret the existing observations. 



In this work, we present a pilot study to explore the properties of hot gas in galactic bulges by 
conducting 3D hydrodynamic simulations. In these simulations, SNe are randomly generated that 
statistically and spatially follow the stellar light distribution. Based on the mean temperature and 
density of the surrounding medium, we adaptively determine the appropriate sizes of individual 
SNRs and generate their density, temperature, and velocity profiles from a library of ID simulated 



^Hereafter in our paper we use the term wind specifically for a supersonic outfiow and outQow for an outflow in 
either supersonic or subsonic state. 
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SNR templates. We then plant such structured SNR seeds into the 3D simulation grid and let 
them evolve. We terminate the simulation when it reaches a statistically steady state. The 3D 
simulations not only provide us with the dynamical structures of SNRs, but also enable us to trace 
the thermal and chemical states of the bulge wind material. 

The organization of the paper is as follows. In §2 we describe the main physical ingredients of 
the bulge wind model and the numerical methods. The results are presented in §3 and discussed 
in §4. We summarize our results and conclusions in §5. 



2. Model and Method 
2.1. Model Basics 

We model hot gas inside a galactic bulge that originates from continuous stellar injection in 
the form of stellar winds, and is heated by sporadic Type la SNe. The dynamics of the hot gas is 
described by the following equations: 



^ +2V ■ {pv) = p,{r) + psM{r,t), 



dpv 
dpE 

P = nkT 



■f V • (pvv) + VP = -pV$, 

+ V • [{pE + P)v] = -pv • V$ + 5s^(r, t) + S,{r) - ntneA{T), 



(1) 

(2) 

(3) 
(4) 



where p, v, P, T, and E denote density, velocity vector, pressure, temperature, and total specific 
energy of the hot gas; and ^ is the gravitational potential field; nt and Ug are the number density of 
ions and el ectrons; n = rit + np^ is t he to tal number density; A(r) is the normalized cooling function 



taken from 



Sutherland &: Dopital (|l993l ). assuming an optically thin plasma with solar abundance; 



and denote the mass and energy input from evolved stars. The mass and energy input from 
individual SNe, /Os]v(r,i) and 5's]v(r,i), are explicitly expressed as a function of position and time. 

We adopt parameters appropriate to the bulge of the Milky Way. The stellar mass d istribution 
of the bulge follows the potential-density pair of the Hernquist profile (jHernquistlll990l ): 



^bulgel,' ) — ; ) Pbulge{' ) 

r + rs 



27r r (r + r^)^ ' 



(5) 
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where is the scale radius and M^^i^^ is the total n i ass of the bulge. 



r, to be 2.4 x 10^" M ^;, and 0.42 kpc (e.g., 



Lepine &: Lerov 



200 



3 



Kent 



1992 



Zhao 



199 



3ere we set the M 



Blum 



1995 



and 



Wolfire et al. 



1995 



Other components of our Galaxy, such as the disk and dark matter halo. 



have only little effects on the gas dynamics within the Galactic bulge and are thus ignored in our 
simulations. 

The stellar mass loss from evolved stars, p*{r), follows the stellar mass distribution. The total 
stellar mass loss rate of the bulge (M = J pdV) is constrained by current theoretical predictions 
and related observations, although the true value cannot be observed directly inside the bulge. It is 
inferred from our knowledge of the stellar population and evolution, together with observations of 
the mass loss of similar stars in the solar neighborhood. Estimates of the stellar mass loss rate may 
vary by more than a factor of two . Assu ming a single stellar population with a standard Salpeter 
initial mass function, ICiotti et al. I (jl99lh found that the stellar mass loss rate can be approximated 



as 



M = 0.25Liotro'-^MQ yr-i, 

where Lio is the current optical blue-band luminosity in units of IO^^L^^q of the stellar 
and tiQ is its age in units of 10 Gyr. Adopting a blue-band luminosity of 2 x lO^L^ . ^:^ i 



(6) 



popu 


ation 


Cox 


2000. 



p571) and an age of 10 Gyr for the bulge, we have M ~ 0.05 Mq yr'^. iMarastonI (20051. fig. 22) 
directly relates mass loss to the total mass of a stellar population (see also lTang et al. II2008I ). which 
gives M ~ O.OTMQyr""*^. Another estimate based on observations of asymptotic giant branch stars 



gives M = 0.64Lio MQyr ^, consiste nt with the ste . 



ellipticals from mid-IR observations (jAthey et al. 



lar in ass loss rate inferred for a sample of nine 



2002). Thus, M for the Galactic bulge can be 



as high as 0.13 yr ^. We therefore run models with different mass loss rates, as discussed in 
52.3. 



The energy feedback of the stellar bulge is dominated by the mechanical energy from Type 
la SNe. Following the convention, we assume that each SN releases 10^^ ergs mechanic energy. 
We take the Type la SN rate of E-SO galaxies to be 0.1 2 SNu (SNu is defined as one SN per 



10 Lb,0 per century; ICappellaro et al. 



19991 : 



Cox 



2000l . p467), which gives about one SN per 



3000 year in our Galactic bulge. The energy input from the stellar wind, 5*(r), is assumed to be 
thermalized to their stellar kinematic temperat ure T^, = unipa'^ /3k ~ 3 x 10^ K, corresponding to 
a stellar velocity dispersion around 100 kms~^ ( Eckart et al. Ill993l ). Overall this energy is almost 
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negligible compared with that from SNe. 



Each SN also ejects an adopted (Chandrasekhar) mass of 1.4 Mq. Though the total amount 
of SN ejecta is much less than that of the mass loss from evolved stars, the SN ejecta contribute 
most of the metals, especially iron. In order to trace how these iron-rich ejecta mix with the stellar 
wind material, we additionally incorporate a separate advection equation 

dpx 



dt 



+ V • (pxiv) = 0, 



(7) 



where Xi is the mass fraction of the ith component, with the constraint = 1- I^i the present 

simulations we have two component s, the iron mass and the rest of gas mass. The iron mass from 
SNe, assumed to be 0.7 per SN (INomoto et al. Ill984l ). is part of the SN ejecta. The iron mass 
fraction in stellar wind material is /fb,© = 0.1% (i.e., the nominal solar). The iron abundance can 
trace the SN enrichment. Any zone with iron abundance greater than the solar value is enriched 
by SNe. Thereafter we refer the iron ejecta as purely from SNe. 



The simulations are performed with FLASH (jFrvxell et al. Il2000l ). an Eulerian astrophysical 
hydrodynamics code with the adaptive mesh refinement (AMR) capability, developed by the FLASH 
Center at the University of Chicago. FLASH solves the Euler equations and uses the piecewise- 
parabolic method to deal with compressible flows with shocks. We take the advantage of the AMR 
capability to accurately include the heating of individual SNRs. 



2.2. One Dimensional Model 

To help set up the 3D simulations, we first simulate a ID model. Assuming that the energy and 
mass inputs are continuous in time and spherically symmetric in space, we may simplify Eqs. (1)- 
(3) into a ID problem. For a specific galactic bulge, the crossing time — the time required for the 
gas to flow from the center to the outer boundary — is a few million years (assuming a size of a few 
kpc for the bulge), significantly shorter than the evolutionary time scale of the stellar energy and 
mass input rates. Therefore, the problem can be regarded as time-independent. The bulge outflow 
can reach a steady state if radiative cooling does not affect the dynamics. 

We run the ID simulations by using a gas-free initial condition and by continuing to add 
corresponding energy and mass in each zone. FLASH handles the mass and energy inputs with 
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an operator split method. At each time step it first solves the Eulerian equations without the 
source terms, then explicitly updates the solution to account for the corresponding source terms. 
To properly conserve the mass, momentum, and energy, we implement this procedure in three 
steps: first, we update the gas density according to the amount of mass input in that step; next, we 
modify the gas velocity to satisfy momentum conservation; finally, we modify the gas temperature to 
conserve the total energy. We verified that this im plementation can exactly produce the analytical 



solution of a star cluster wind ( Canto et al 



hood ). 



The system eventually evolves to a stea dy state. Such a steady o utflow solution can be an- 
alytically derived without including cooling (jWhite &: Chevalier] Il983l ). If the ID bulge outflow 
solution has a sonic point, a subsonic outflow can then develop into a supersonic outflow (i.e., a 
bulge wind). The final state of such a wind does not depend on the specific initial condition. Under 
certain conditions (e.g., due to significant radiative cooling or low specific energy input; see §4.2 
for more discussion), the sonic point may not exist and the gas outfiow may be sensitive to the 
boundary condition as well as the initial condition. The use of the outfiow (sometimes called zero- 
gradient) boundary condition would then introduce an artificial force inserted by the leveled-off 
pressure that would produce perturbation, propagating inwards on a time scale comparable to the 
crossing time. This situation would also occur if a simulation region were too small to include the 
sonic point (if present). Thus we use the ID solution to make sure that the sonic point is included 
in the 3D simulation domain (a cubic box). 



2.3. Three Dimensional simulations 

Two 3D simulations are performed to examine the properties of galactic bulge winds. Their 
key parameters are listed in Table [1] for quick reference. The major difference between these 
two simulations is the mass loss rate: 0.05 Mq yr~^ for Model A (the 3D reference model) and 
0.1 Mq yr^^ for Model B, representing the uncertainty in the mass loss rate (§2.1) or extra mass 
loading expected in galactic bulges (§4.2; see also Li & Wang 2009 in preparation). The highest 
spatial resolutions are ~ 3.9 and 4.9 pc respectively for the two models. The effective single-grid 
resolution of the simulations are 1024^ and 2048^ zones. The steady wind fiow established in the 
ID model is used as the initial fiow with an iron mass fraction of the solar value (/ = 0.1%). 
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In the 3D realizations, SNe explode randomly according to a Poisson process with the mean 
overall rate. Their spatial distribution statistically follows the stellar density distribution. It would 
be computationally very expensive, if even possible, to simulate the evolution of each SNR on 
sub-parsec scales within the bulge-wide flow. Instead, we adaptively plant individual structured 
SNR seeds into the 3D simulation grid and then let them evolve. We do not simply adopt the 
Sedov solution, which is generally not appropriate for an SNR evolving in a hot tenuous medium 



(|Tang &: Wangi |2005| ) . especially when the dynam ics of the SN e . 



to a scaling scheme detailed in a separate paper (jTang &: Wang 



ecta are considered. According 



2009|), the structure of an SNR 



can be scaled from a template SNR simulated in a different ambient medium setting. We have 
constructed a library of template SNRs from ID simulations, assuming a selection of ambient gas 
temperatures and densities. Each entry of this library consists of the density, temperature, and 
velocity profiles at a particular age and a forward shock radius. 



We apply the scaling scheme to dynamically generate the profiles of each SNR seed. Specifically, 
we select a spherical region around each SN location within whic h the density and pressure are 



sufficiently smooth, using the Loner's error (FLASH User's Guide; iLoner 



19871 ) as the estimator. 



We find that at least 500 zones (i.e., more than 5 points for the radial profiles) are needed to 
reasonably well represent a structured SNR seed. This in turn requires the minimum embedding 
radius to be at least 20 pc, given the spatial resolution that is achievable in our simulations. With the 
embedding radius determined, we calculate the mean density and gas-mass-weighted temperature 
of th e enclosed gas to find the most suitable template in the library and to form the required SNR 



seed (ITang k, Wang 



200 



1 



The planting of an SNR seed also takes a few steps. First we refine the affected region to 
the highest refinement level. Then we normalize the SNR structures to ensure the conservations 
of mass, momentum, and energy within that region (see Appendix A for details). Finally, the 
innermost region that encloses 0.7 Mq is traced as the pure iron ejecta of the embedded SNR. This 
embedding procedure is well parallelized and allows for linear scaling up to at least 1024 processors. 



To save computing time, one may simulate only one octant of the bulge by adopting a refiecting 
boundary condition at the surfaces across the center. A test run, however, shows that a refiecting 
boundary condition introduces correlated wave interactions when an SN explodes near the refiecting 
boundaries. This effect is not physical and difficult to quantify. It is most serious near the bulge 
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center where three reflecting boundaries intersect and the stellar density, hence the rate of the SNe, 
is the highest. Thus we resort to simulating the whole bulge, which is centered inside the simulation 
domain. However, we only simulate one octant at full resolution, while the highest resolution in 
the rest of the grid is degraded by a factor of four, except for regions where SNRs seeds have just 
been embedded. These regions are forced to have the full resolution in all octants, and are held for 
10^ years before returning to the default AMR. We use refinement estimators acting on the density 
and pressure to determine whether a block needs to be refined or derefined, adopting the default 
criterion suggested (in FLASH User's Guide, i.e., refining a block if any estimator is greater than 
0.8 and derefining it if all estimators are less than 0.2). Regions outside the sonic radius (which is 
obtained from the ID model) are allowed to have a lower refinement level that gradually decreases 
with radius. This approach circumvents the reflection boundary problem at the expense of about 
60% more computing time, which is acceptable. In addition, it allows us to examine the resolution 
effect within a single run. 

A statistically steady state of such a 3D simulation can be reached after a few crossing times. 
We quantify the establishment of the steady state of a 3D bulge wind by examining the relative 
variation of its global quantities such as the total mass and energy. Here we define the variation 
as the change of a given quantity relative to its initial value, i.e., the relative difference between 
3D and ID. The variation of the total mass within 2.0 kpc of Model A is shown in Fig. [T^ by the 
solid line. Compared to its initial value, the total mass increases to ~ 7% on average and fluctuates 
around this value with a period of '^lOMyr, comparable to the flow crossing time. As expected, 
the mass variation within the inner 1.2 kpc radius, displayed as the dashed line in the figure, has 
a shorter fluctuation period. The expected iron mass fraction, if fully mixed with stellar wind 
material, is 0.35% 3.5 times the solar abundance) in Model A. We show the variation of the 



Table 1. Model Parameters 



Model 


M 


Esn 


^sonic 


AL 


L 




(Mo yr-i) 


(10'*°ergss"') 


(kpc) 


(pc) 


(kpc) 


A 


0.05 


1.1 


1.0 


3.9 


4 


B 


0.1 


1.1 


1.8 


4.9 


10 
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iron abundance (in the solar units) in Fig. \T}p. It takes about 5Myr for the hot gas inside 2kpc 
to gain the expected iron mass. The variation of the iron abundance is smaller than that of total 
mass, mainly because it actually only reflects the ratio of total iron mass to the total gas mass. 
By introducing random SN events in the 3D simulations, the globally conserved quantities are no 
longer constant as they should be in a ID spherical steady flow. Only when a hydrodynamic steady 
state is established in the simulations (i.e., the fluctuation has reached a statistically stable level) 
is the comparison between ID and 3D results meaningful. 

We check the resolution effect based primarily on X-ray luminosity, which is particularly sen- 
sitive to the density structure in the simulations. We find that the X-ray luminosity difference 
between the high resolution octant and the other seven octants is rather small. This is partly 
because the majority of the emission arises from individual SNRs which are resolved at the same 
resolution in all the octants. To examine the resolution effect more directly, we resume Model A 
with an increased spatial resolution by a factor of two, and let it only evolve for 0.1 Myr, limited 
by the available computing time. This simulation produces finer structures of the bulge gas, and 
the resultant X-ray luminosity increases about 3% in the high resolution octant and about 10% in 
the other seven octants. This demonstrates that our results are quite robust and are only slightly 
affected by the spatial resolution. 



3. Results 

In this section we present the gas properties extracted from the 3D hydro dynamical simula- 
tions. We first detail the results of Model A (Fig. [2|) and then present Model B (Fig. [3]) for com- 
parison. Data near the outer region are excluded in our analysis to avoid any potential artifacts 
introduced by the assumed outer boundary condition of the simulations. We show time-dependent 
gas properties such as global structures, individual SNRs, and X-ray luminosities as well as various 
time-averaged measurements. The average is made over a time span ranging from 15 to 30 Myr, 
when the simulations have reached quasi steady states (see Fig. [1]). 
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3.1. Structures 



Fig. [5] shows snapshots of the simulated density, temperature, pressure, and iron ejecta mass 
fraction of Model A in the z = 2 pc plane. Sporadic SN events produce non-uniformity in the bulge 
wind. The prominent features are various shell-like and filamentary density structures. These shell- 
like structures of SNRs are easily identified in the outer region, where the explosions are less frequent 
and each SNR can evolve individually to a large volume before colliding with others. Individual 
SNRs near the bulge center appear more compact because of the high gas density and pressure 
and frequent interactions with adjacent remnants. Evolved remnants tend to be dispersed and 
advected outward. Higher temperature regions in general represent low-density interiors of SNRs. 
The distribution of the pressure is much smoother than those of th e density and temperature 
(Fig. [2li) , as has also been f ound for the SN driven IS M in the Galaxy (Ide Avillez &: Breitschwerdt 



2005 



Mac Low et al. 



2005 



Joung fc Mac Lowl 120061 ) and in starbursts ( 



Joung et al 



20081 ). The 



shell-like structures in the pressure map correspon d to the expan ding blastwaves. Inside each SNR, 
the pressure is nearly uniform, as expected from the lsedoy" ( 1959 ) solution. The spatial distribution 
of iron ejecta is far from uniform, as illustrated in Fig. (2}:. Regions with the lowest iron mass fraction 
are primarily filled with stellar wind material with little mixing with the iron ejecta. Gas with an 
intermediate iron mass fraction represents iron ejecta diluted by constantly injected stellar wind 
material (or mildly by numerical mixing). Though being diluted and fragmented, the iron ejecta is 
advected out of the bulge, hardly mixing with the bulk of stellar wind material. Hence, the ISM is 
not uniformly enriched by SN ejecta within the bulge. Similar results are also present in Model B 
(see Fig. El). 



Fig. [3] shows sample density, temperature, and velocity profiles of the hot gas at two rep- 
resentative times. Individual troughs shown in the density profiles represent interiors which are 
surrounded by peaks (i.e., expanding shells of SNRs). The temperature profiles are nearly anti- 
correlated with the density profiles: a peak in temperature usually corresponds to a trough in 



density at the same locus, whi ch is equivalent to the smooth pressure profiles (e.g.. 



Mac Low et al 



2005 



Joung fc Mac Lowl 120061 ). As SNRs evolve, the loci of peaks and troughs change with time. 



Multiple evolved SNRs appear to be wave-like. This wave-like flow driven by sporadic SNe produces 
the fluctuations in the globally conserved quantities (Fig. [T|) . 



We demonstrate the evolution of a few SNRs (labeled as I, E, and IE) very close to the bulge 
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center in Fig. [5] taken from the high resolution study. The forward blastwave of SNR I is evident 
in the density panel at 10"^ year. The pure iron core has a radius of about 10 pc at this time. 
The blastwave expands faster toward the lower-right region, where a lower density cavity has been 
created by an earlier SNR (labeled as I). SNR I is unusual in that it occurs right at the bulge 
center. Because of the high stellar wind injection rate there, its iron ejecta are diluted quickly. 
SNRs away from the center are less affected by the stellar wind dilution. SNR E, for example, 
at an age of 2 x 10^ year still has an iron mass fraction of ~ 6%, corresponding to 60 times the 
solar abundance. The collision between SNR I and II is also evident in the lower-right subpanel of 
each group. The evolved shapes of these SNRs are asymmetric because of both the inhomogeneous 
environment and interactions with other SNRs. 



3.2. Time- Averaged Distributions 

Fig. [6] shows the time-averaged gas mass and volume distributions in the high resolution octant 
as functions of temperature and density in two regions: < r < 0.6 kpc and 0.6 < r < 1.2 kpc. 
The distributions in both regions are broad. As expected, the mass distribution is biased toward 
the lower temperature and higher density side relative to the volume distribution. Relatively cold 
dense gas (e.g., distributed in the upper-left corner of Fig. [6^) occupies negligible volume, lying 
outside the 1% contour level of the gas volume distribution. To see the results more quantitatively, 
we show their marginal distributions in temperature and density respectively (Fig. [7]). The volume 
distribution resembles a power law with an index of roughly —2 at the high-temperature end. The 
mass distribution is similar to the volume distribution but peaks at a lower temperature. 

Fig.Oshows the differential emission measure (EM) as a function of temperature within a radius 
of 1.2 kpc. The EMs from the low- and high-resolution regions are nearly identical around the peak 
value. The bulk of the broad EM distribution, peaked at ~ 3.7 x 10^ K, can be approximated by 
a log-normal distribution, with small deviations mostly at the high temperature end. The grey 
region in the figure encloses the 50% intervals below and above the mean of the whole region. 
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3.3. X-ray Emission and Spectra 

We calculate the luminosities and spectra using the standard software package XSPEC based 
on the EM distributions. For the hot gas considered here, we adopt the standard MEKAL model 



( Mewe et al 



1985 



Liedahl et al 



19951 ) for plasma in collisional ionization equilibrium. Fig. shows 
a synthesized spectrum of Model A assuming a solar abundance. The luminosities in a few bands of 
Model A are listed in the 4th row of Table [21 The corresponding spectra and luminosities of Model 
B are presented in Fig. [9] and in Table [2] as well. The luminosity of Model B in the low energy band 
increases dramatically; e.g, in the 0.2-0.5 keV band it is more than 60 times larger than that of 
Model A. This increase is due to the combination of higher density and lower temperature (hence 
higher emissivity) of the weak bulge wind (see §4.2 discussion). Of course, the X-ray emission 
depends on the metal abundance as well. If an abundance expected for stellar material fully mixed 
with SN ejecta were adopted, the luminosities in the three bands would increase by a factor of a 
few, as listed in Table [2j 

Since the iron ejecta are in fact not well mixed with the surrounding material, we examine the 
effect of the non-uniform metal distribution on X-ray emission. The normalized EM distribution 
as a function of the iron abundance is plotted in Fig. [10] as the solid black line. It shows that 
about 80% of the EM comes from the material with the iron abundance less than 1.5 times solar. 
The corresponding distribution of the luminosity in the 0.3-2 keV band nearly follows that of the 
EM. But the distribution of the luminosity in the 2.0-5.0 keV band is affected considerably by 
gas with higher iron mass fractions (corresponding to SNR interiors that in general have higher 
temperatures). Thus the majority of the X-ray emission from the bulge wind comes from the 
stellar wind material that is hardly enriched by SNe ejecta. In the following, we thus adopt the 
solar abundance for our calculations. 



Fig. [H] shows the time variation of the 0.3-2.0 keV luminosity of Model A in two concentric 
regions: an inner region with r < 0.6 kpc, and an outer region with 0.6 < r < 1.2 kpc. The 
luminosity in the inner region has a larger fluctuation (up to a factor of 3) than that in the outer 
region. The X-ray emission from the region with r > 1.2 kpc is negligible. The total 0.3-2.0 keV 
X-ray luminosity is only ~ 10^^ ergss^"*^ with a fluctuation less than a factor of two. Thus overall 
the X-ray emission of the diffuse gas is not significantly affected by the sporadic SN heating. 
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Fig. [12] illustrates the 0.3-2.0 keV surface brightness profile at several representative times. 
The brightness varies by more than an order of magnitude near the center. Compared with the 
stellar surface density profile (gray line), the X-ray profiles are generally steeper. Fig. 1131 shows 
the X-ray surface brightness maps in three representative bands (0.3-0.7, 0.7-2.0, and 2.0-5.0 keV). 
The maps in the two lower energy bands do not show significant structures; those small features 
are typically not associated with individual SNRs. In the 2.0-5.0 keV map, individual SNRs are 
recognizable because they are the primary source of the hard X-ray emission. 

3.4. Energy Distribution 

From the 3-D realization, we can further quantify the thermal, kinetic energies as well as the 
turbulent motion of the hot gas. Fig. [14] shows the distributions of kinetic and thermal energies 
as a function of temperature at four different regions. The energy distribution peaks at 5 x 10^ K. 
At temperature below 5 x 10^ K, the energy is dominated by the thermal energy near the center 
(0 < r < 0.4 kpc), while the thermal and kinetic energies contribute almost equally at large radii 
(r > 1.2 kpc). Note that for ideal gas with the ratio of specific heat 7 = 5/3, the Mach number 
is 1.34 (i.e., supersonic flow) when its thermal energy equals to its kinetic energy. The thermal 
energy of much hotter gas, primarily low-density hot SNR ejecta, is always much larger than its 
kinetic energy, which means that the SNR ejecta leave the simulation region subsonically. The 
difference between the total kinetic energy and its radial component reveals the energy contained 
in non-radial motion. Inside the scale radius (0 < r < 0.4), ~ 30% of the kinetic energy is in 
non-radial motion. This fraction is less than 2% in the outer shell (1.2 < r < 1.6). As a whole, > 
80% of the thermal and kinetic energy of the bulge wind is stored in gas with temperature between 
10^-^ K to lO'^-^ K. The hotter gas (T > lO'^-^ K) contains < 5% of the total energy. 

4. Discussions 

4.1. Comparison between the ID Model and 3D Simulation 

Fig. [15] compares averaged radial profiles for the density, temperature, velocity, and pressure in 
a few snapshots of Model A with the corresponding ID results. The velocity and pressure profiles 
averaged in radial bins follow the ID results excellently. In the 3D simulation there is no unique 
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location of the sonic point to divide the supersonic flow from subsonic flow since the velocities 
and temperatures vary greatly from point to point. Nevertheless, the average hot gas outflow does 
become supersonic at ~ 1.0 kpc, close to the ID sonic point (marked by the arrow). The density 
and temperature profiles deviate significantly from the corresponding ID results in the inner region 
(r < rg). In particular, the density profile tends to be more centrally peaked in the 3D simulations. 
This temporary accumulation of stellar wind material is largely due to the continuous mass injection 
and the lack of prompt heating. The variation of the density profiles, especially near the center, 
is closely coupled with the realization of SN events. At large radii the density and temperature 
profiles become nearly identical to the ID results. Since the radiation of the bulge wind is negligible, 
the identical radial profiles at large radii are expected due to the conservation of mass and energy. 
Thus the nature of sporadic SN explosions does not affect the overall gas dynamics on large scales. 
In Fig. [16] we compare the profiles of the octant at full resolution to profiles of three low resolution 
octants. The figure suggests that the profiles show numerical convergence at a level better than 
the random variance among realizations. 

On average the gas temperature in the 3D simulation has a lower and flatter profile than that 
in the corresponding ID result. In the ID model, the gas temperature is the highest at the center 
and decreases outwards monotonically. In contrast, the average gas temperature in a 3D simulation 
can be much lower at the center than in the surroundings, because of non-uniformly distributed SN 
heating. The distribution of SN heating depends on the ambient medium. An SN that explodes in 
a dense environment like the bulge center, heats a small region to a rather high temperature. The 
small amount of overheated gas advects outward and carries a large fraction of the SN energy with 
it, while leaving much of the central gas unheated. The unheated gas accumulates around the bulge 
center, resulting in a relatively steep density profile near the bulge center (see Fig. 1151 panel a). The 
large density gradient makes it even harder for the SN heating to be unif ormly distributed near th e 



center, so it tends to be transported to outer low density regions (see also lHnatyk &: Petruk Ill999l ) 



Thus a low-temperature inner region naturally forms under the sporadic SN heating scenario. 

It is also worth noting that the temperature profile depends on the weighting methods, as shown 
in Fig. [TTI The EM-weighted value, which closely resembles that inferred from X-ray observations, 
is relatively low because it is primarily determined by the low-temperature denser material (see 
§3.2). In comparison, the temperature profiles from different weighting methods are identical in 
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the ID model. This is because in the ID model, all quantities such as temperature, density, and 
pressure are monotonic functions of radius. Therefore the density has a one-to-one correspondence 
to the temperature, which also explains the single-line gas distribution of the ID model in the 
temperature and density space (Fig. [6]) . 

The integrated gas properties of Model A are also significantly different from those of the ID 
model, such as the EM distributions, spectra, X-ray luminosities, etc. In the ID model the EM peaks 
at ~ 8 X 10^ K and is truncated sharply at low or high temperatures (three-dots-dashed magenta line 
in Fig. [8]), while the EM distribution of Model A peaks at a lower temperature and is much broader. 
But the difference between the spectral shape of Model A and the corresponding ID model is small 
in the 0.7-3.0 keV band (see Fig.[9|). Model A has slightly higher (about 30% more) X-ray luminosity 
in this band. At lower (< 0.7 keV) and higher energy bands (> 3.0 keV), Model A gives much higher 
photon fluxes, especially for some line features, such as OVII (22. oA, 0.56 keV) and helium-like iron 
Ka (FeXXVKa; 6.7 keV), largely due to the much broader temperature distribution. It is worth 
noting that the FeXXV Kq line emission in Model B increases but the line in its corresponding ID 
model is missing because of the relatively low and narrow gas temperature distribution in the ID 
model. 



4.2. Observational Implications 

There are two direct predictions from the 3D simulations that help to understand partly the 
observational puzzles faced by the ID wind model. First, the relatively low emission-weighted gas 
temperature is consistent w ith the results inferred from the diffuse X-ray emission in g alactic bulges 



and elliptical galaxies (e.g. 



Sarazin et al 



2001 



David et al 



20061 : 



Li fc Wangll2007l ). Second, the 



SN la ejecta are not fully mixed with the stellar wind material inside the galactic bulge, and most 
of the X-ray emission is contributed by the stellar wind material shocked by SN blastwaves. Thus 
the X-ray spectra in general reflect only the metal abundance of the stellar wind material, which 
is consistent with the app arent solar metallicity infe rred from the diffuse X-ray emission in a large 



sample of ellipticals (e.g.. 



Humphrey &: Buote 



20061). 



In addition, the broadening of the gas temperature distribution can also affect the determina- 
tion of metallicity. A low metal abundance could be obtained by fitting the X-ray spectra of such 
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gas with an isothermal plasma model. This effect is demonstrated in Fig. [THJ We simulate the 
X-ray spectrum based on the MEKAL model in XSPEC, adopting the approximate log-normal EM 
distribution (see §3.2) and a solar abundance. The data points in the upper panel show a simu- 
lated spectrum of about 600 photons. If we fit this simulated spectrum with a single-temperature 
MEKAL model, the resulting abundance is only about half solar (0.49 it 0.17). The fit is statisti- 
cally acceptable (with ^ = 42.4/41), partly due to the small counting statistics. If the simulated 
spectrum contains 10^ photons instead, the fit is no longer acceptable (with = 213/91) and 
the best fit tends to give a much lower abundance (about one-third solar). A two-temperature 
component model can significantly improve the fit, and the abundances are in general less than one 
solar, although they are not strongly constrained. Thus the abundance is likely underestimated 
by fitting an X-ray spectr um of gas with broad temperature and density distributions (see also 



Strickland &: Stevens 



19981) 



Given the low resolution and small signal-to-noise ratio of X-ray spectra typically available 
for diffuse hot gas, it is difficult to distinguish gas with a broad temperature distribution from an 
isothermal plasma, especially within a narrow energy band. In Fig. [9]we plot an arbitrarily normal- 
ized spectrum of a single-temperature hot plasma (0.8 keV for Model A and 0.4 keV for Model B, 
where the EM distribution of the corresponding ID model peaks). The spectra of the 3D simulation 
and the isothermal plasma model are similar in the 0.5-2.0 keV band. This approximation gives 
a potential shortcut to fit the X-ray observations coarsely with only one or two gas components, 
which are not necessarily able to reveal the actual physical and chemical state of the ISM. With 
higher resolution spectra, line diagnostics may provide additional useful information to reveal the 
gas properties. 



However, the models still predict far lower X-ray luminosities than observed. The observed 
diffuse X-ray l uminosity in the . 5-2 keV band of our Galactic bulge and M31 bulge is abou t 



XO^s ergs s ^ ( Shirev et al 



2001 



Takahashi et al. 



20041 : 



Muno et al. 



20041 : iLi fc Wang I l2007l l 



Using the reference mass and energy input rates, the bolometric X-ray luminosity predicated by 
the ID wind model is no more than 10'^^ ergs s~^. The inhomogeneous structures of the gas density 
and temperature in Model A increase by no more than half an order of magnitude the luminosity in 
the 0.5-2 keV band (Table 2). This under-prediction of the observed luminosity remains a serious 
problem for the models. 
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The parameter with the strongest influence on the luminosity of the bulge wind is the stellar 
mass loss rate. As demonstrated in Model B which has a doubled mass input rate, its mean gas 
density increases while the temperature decreases. The X-ray emission increases by a factor of ~ 20 
in the 0.5-2.0 keV band. This enhancement is due to the increase of both density and emissivity. 
To quantitatively understand this trend, let us first consider a scaling relation of the bulge wind. 
The specific heating (3 = E/M determines the mean central gas temper ature of the wind solution. 



In a steady flow, the velocity u oc (3^'^ (e.g. IWhite &: Chevalier! Il983l ). so long as gravitational 



potential energy remains unimportant. We have p oc M[5'~^'^ from the mass conservation equation, 
and T (X (3 from the energy conservation equation. For gas with a temperature between 0.5 
and 1.0 keV, the X-ray emissivity in the 0.5-2.0 keV band is roughly inversely proportional to the 
temperature [e.g., A(r) (x T~^''^ based on the MEKAL model]. The total X-ray emission is then 
L oc p'^A{T) oc M^'''£^~^'''. For winds with the same /? the velocity and temperature proflles are 
the same, and density proflles differ only in their normalization, which is proportional to M. We 
thus use to denote the separate change of either M or E, e.g, L oc M^''' is equivalent to 

L oc (/?//?* )~^'^ for a flxed energy input, and L oc E~^-'^ to L oc (Z?//?*)"^'"^ for a flxed mass input. 
Fig. [T9l shows the luminosity in the 0.5-2 keV as a function of (3 based on a suite of ID simulations. 
In this plot the luminosities are normalized to that of a ID reference model (L* ~ lO'^^ergss"^; 
corresponding to Model A) . The scaling relations of L versus M and E closely match the simulations 
when is larger than 0.5. In the adopted galactic bulge, if the speciflc energy approaches one- 
third of the reference value, the simulated luminosities increase sharply and no longer follow the 
scaling relation, because the gravitational potential becomes dynamically important for such a small 
f3. The corresponding results of Models A and B are also plotted in Fig. [19] to show the effects of 
gas inhomogeneity. The ratio of the 0.5-2.0 keV luminosity of Model A to that of the corresponding 
ID model is about 2, and the ratio is about 3 for Model B. 

Let (3 be flxed at the reference value in Model A, the velocity and temperature proflles of the 
wind should remain the same and Lx oc Af^. Hence, both E and M need to increase by an order 
of magnitude in order to boost the X-ray luminosity to match that observed in the M31 bulge. 
Although a bulge wind with a reduced (3 can signiflcantly increase the diffuse X-ray emission as 
demonstrated in Model B, it is unlikely to be the best solution. In the presence of a reasonable 
dark matter halo with properties like that of the Milky Way galaxy (e.g., M^aio — IO^^Mq, Vmriai ~ 
250 kpc, an NFW proflle with a concentration parameter of 15), ~ 120% of the available energy 



- 19 - 



input in Model B is required for the escape of all the stellar ejecta to beyond the virial radius. 
Thus the bulge wind in Model B is not energetic enough to escape the galaxy potential well and 
the bulge wind material must accumulate within the virial radius. This accumulated material may 
form a considerabl e circum-galactic medium (CGM) that could eventually quench the bulge wind 



(jTang et al. 



200. 



The gas outflow might be subsonic in the vicinity of galactic bulges and other ellipticals with 
low Lx/Lb ratios under the interaction between the bulge wind and CGM. The subsonic state 
is necessary to explain both the large luminosity (compared to the wind model prediction) and 
the extent of the diffuse X-ray emissions. Even if a bulge wind has the energy to escape the halo 
virial radius, the wind can be reversely shocked and stalled by the CGM. When the reverse shock 
propagates inward within the sonic point, the bulge wind turns into a globally subsonic outflow. 
Such a subson ic state can be qua si stable for a long time with proper treatment of the boundary the 



galactic flow (jTang et al. 



20081 ). A proper 3-D simulation including the entire CGM is possible, 



though computationally very expensive at present. 

Other effects, which are ignored here for simplicity, can further affect the X-ray luminosity of 
the galactic bulges. A vertically collimated bulge wind can be a promising way to significantly boost 
the X-ray luminosity, particularly for objects such as M31 with a considerable mass in the galactic 



disk. Motivated by the observed bipolar diffuse X-ray emission in the M31 bulge (jLi &: Wang 



20071 ). we have found in preliminary simulations that a vertically collimated wind, confined by 
the surrounding gaseous disk, will have significantly higher X-ray emission (by up to an order of 
magnitude). The major reason is that the density diverges much slower in the collimated wind 
than in the spherical wind. We plan to address this issue in a separate paper (Tang & Wang in 
preparation). 



5. Summary 

In this work, we have explored the properties of the structured hot gas created by sporadic 
SN explosions inside a galactic bulge by conducting detailed 3D hydrodynamical simulations. Our 
main results are as follows: 
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• A galactic wind may be generated in a galactic bulge with the standard empirical stellar mass 
loss and Type la SN rate. The gas properties fluctuate in time particularly in the central 
region lying within the sonic radius of the wind, where individual SN explosions strongly 
influence the density and temperature distributions. At larger radii, the spherically average 
profiles of the 3D simulations follow those of the ID models. Therefore the ID treatment of 
a galactic bulge wind flow is a reasonable approximation on large scales. 

• Sporadic SN explosions produce 3D filamentary and shell-like structures in the gas. These 
structures result in broad density and temperature distributions, compared to the ID model. 
Furthermore, the relatively low temperature of the structures leads to an emission measure- 
weighted temperature that is significantly lower than the expected value inferred from the 
specific heating and has a relatively flat radial distribution throughout the bulge region, 
consistent with observations. 

• Iron ejected by SNe does not mix well with the surrounding gas within the bulge region and 
has a relatively high temperature and low density, so it contributes primarily to emission in 
the energy band > 2keV. The diffuse soft X-ray emission comes from shells associated with 
SN blastwaves, which are hardly enriched by SN ejecta and have a metallicity close to the 
ambient gas that originates in stellar winds. This, together with the temperature broadening, 
helps to explain the apparent sub-solar abundance of the soft X-ray-emitting gas in galactic 
bulges / ellipticals . 

• Compared with the ID spherical wind model, the structured hot gas in 3D simulations can 
boost the X-ray emission in an intermediate energy band (e.g., 0.5-2.0 keV) by a factor of a 
few. This increase is more significant at the lower or higher energy bands due to the broad 
distributions of both temperature and density. In order to produce the luminosity and surface 
brightness distribution similar to the observed diffuse X-ray emission, the bulge outflow likely 
needs to be in a subsonic state and/or an angularly confined configuration. 
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Time (Myr) 

Fig. 1. — Variation of the total gas mass [panel (a)] and energy [panel (b)] relative to their initial 
values and iron abundance [panel (c)] in Model A. Dash blue lines show the variation within 1.2 kpc; 
other lines show the variation within 2.0 kpc. The three-dots-dash red line in panel (b) denotes the 
variation of the total thermal energy. 
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Fig. 2. — Snapshot of Model A in the z = 2pc plane, showing (a) density, (b) temperature, (c) iron 
mass fraction, and (d) pressure. All plots are logarithmically scaled according to the color bars. 
Note that the upper right quarter region in each panel denotes the data from the octant at full 
resolution. 
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Fig. 3. — Snapshot of Model B in the z = 2.5 pc plane. All plots have the same meaning as in 
Fig. El 
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Table 2. Time-averaged X-ray Luminosities 



Enery band(keV) 


0.2-0.5 


0.5-2.0 


2.0-10 


Iron Ka 


(ergs s""*") 


(lO^^^) 


(10»«) 


(10^^) 


(10^2) 


(A) ID: Zq/Z.^Zq 


0.092/0.22 


0.72/2.42 


0.022/0.066 


0.095/0.31 


(A) 3D: Z0/3.5ZQ 


0.39/1.18 


1.22/4.14 


0.021/0.062 


4.1/13.8 


(B) ID: Z0/I.8Z0 


3.44/5.32 


10.8/18.5 


0.03/0.05 


... / ... 


(B) 3D: Zq/1.8Z0 


26.8/49.8 


29.09/50.4 


0.066/0.10 


7.9/13.5 



Note. — For each model the X-ray emissions are calculated using two 
metallicities: 1.0 and 3.5 solar for Model A; 1.0 and 1.8 solar for Model B. 




-2-10 1 2 

r (kpc) 



Fig. 4. — Two sample realizations of density, temperature, and velocity profiles along x-axis. The 
ID results (gray lines) are included for comparison. 
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Fig. 5. — A sample SNR (centered at I) evolving near the bulge center. The density, temperature, 
iron fraction, and pressure are grouped into four large panels; each panel gives snapshots at four 
ages: 10^, 2.0 X 10^, 4.0 X 10^, 9.0 x 10^ year, ordered from left to right and top to bottom. All 
values are logarithmically scaled according to the color bars. 
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Fig. 6. — Gas mass (solid black) and volume (dash blue) distributions in the temperature-density 
space in two regions: r < 0.6 kpc (left panel) and 0.6 < r < 1.2 kpc (right panel). Contours 
are normalized to their peak values and labelled logarithmically. The shaded regions denote the 
50% coverage for the mass (solid shading region) and volume (hatched shading region) centered 
on their respectively peak values. The corresponding results of the ID model are shown by the 
three-dots-dashed magenta line in each panel. 




6 7 8 -5 -4 -3 -2 

log,o T (K) n (cm-^) 

Fig. 7. — Gas mass (thick solid) and volume (thick dash) distributions as functions of temperature 
(left-panel) and density (right-panel) within r < 1.2 kpc. The corresponding distributions of the 
ID model are shown in thin magenta lines. 
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Fig. 8. — Time-averaged EMs as a function of temperature: the whole simulation region (solid black 
line); the low resolution octants (dash-dot blue line); the high resolution octant (short dash red 
line); and the ID model (three-dots-dashed magenta line). The EMs of the low and high resolution 
regions are linearly scaled to the whole region according to their respective volumes. The grey 
region marks the 50% range around the mean EM of the whole region. Long-dash cyan line is a 
log-normal fitting with the mean log^g ^ = 6.57 and the standard deviation of 0.24. The bin size 
of the temperature is 0.01 on the logarithmic scale. 




10.0 



Fig. 9. — Time-averaged spectra within r < 1.2 kpc (solid black line) of Model A (top panel) and 
Model B (bottom panel), compared with their corresponding ID model (dashed magenta line) and 
arbitrarily normalized spectra of isothermal plasma with a temperature of O.SkeV or 0.4 keV. The 
solar abundance is used to generate all the spectra. 
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Fig. 10. — Distribution of the emission measure (solid line), as well as the luminosities in the 0.3-2.0 
keV (dash red line) and the 2.0-5.0 keV (three-dots-dash blue line) bands as functions of the iron 
mass fraction. 
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Fig. 11. — Illustration of the 0.3-2.0 keV luminosity fluctuation as a function of time in the region 
within r < 0.6 kpc (three-dots-dash blue line), between 0.6 < r < 1.2 kpc (dash red line), and the 
total volume (solid black line). 
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0.7-2.0 keV 2.0-5.0 keV 




-10 1 -10 1 -10 1 

Fig. 13. — Representative intensity maps at full numerical resolution in the three energy bands. 
The unit of x and y axis is kpc. The intensity is in the units of ergss"^ sr~^ and is logarithmically 
scaled. 
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Fig. 14. — Time-averaged thermal (Ei) and kinetic energy (Ei^; E/^^r is the radial component) 
distributions versus the gas temperature in four radial ranges as marked at the bottom of each 
panel. 
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Fig. 15. — Radial profiles of density (a), velocity (b), temperature (c), and pressure (d) sampled 
at four times. The solid lines (f=O.OMyr) show the ID results. The density profile is weighted by 
volume and the other profiles are weighted by mass. The arrow in panel (b) marks the sonic point 
of the ID model. 
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Fig. 16. — Radial profiles of density (a), velocity (b), temperature (c), and pressure (d) sampled at 
four octants. The solid black lines show the results of the octant at full resolution, and other lines 
of three octants of low resolution. 




Fig. 17. — A sample radial temperature distribution weighted by volume (denoted as wv), by mass 
(wm), and by emission measure (wem). The long dash line denotes the ID temperature profile. 



Energy (keV) 

Fig. 18. — Simulated Chandra ACIS spectra of the X-ray emission from the modeled bulge gas, 
based on the differential EM shown in Fig. [8l About 600 and 10^ photons are generated for the 
upper and lower panels, respectively. The solid red line denotes the fit of an isothermal plasma 
model while the dash-dot blue line denotes the simulated spectrum. The normalization is arbitrary. 
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Fig. 19. — Dependence of the 0.5-2.0 keV luminosity on (3. The Hnked diamonds denote the results 
with a fixed M, while the crosses mark the results with a fixed E. The 3D simulations results are 
marked by the star for Model A and the triangle for Model B. Gray lines denote the corresponding 
scaling relationship (see §4.3). 
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A. SNR Embedding Scheme 



After all quantities of zones within the affected region are replaced by the dynamically gener- 
ated SNR profiles, we make additional adjustments to conserve the mass, momentum, and energy, 
accounting for the geometrical effect due to discrete zones and the possible net bulk motion of gas in 
that region. Some symbols are introduced in the following adjustments for clarity: "M" denotes the 
total mass; "IvlVX", "MVY" and "MVZ" are momentum in X, Y and Z directions, respectively; "EK" 
is the total kinetic energy; "EI" is the total thermal energy; the subscripts "0" and "1" denote the 
values before and after the SNR seed is planted, "2" stands for the adjusted values. The symbols 
followed by (i,j,k) denote the values at each zone. In all the calculation a sub-zone method (e.g., 
subdividing a zone to determine whether a zone is partially within the affected region) is used in 
order to achieve high accuracy. 

(a) Conservation of mass. It is accomplished by simply multiplying a factor, fdens = M2/M1 = 
(Mq + Mejccta)/Mi, to the density of scaled SN profile: 

Piihj^k) = fdens X Pi(i,j,A:) . (Al) 
Mo = J^Po(^i,^)^^M> (A2) 
Mi = Y,pAhj.k)^Vi,j,k (A3) 

Mejecta denotes the SN la ejecta and is 1.4 Mq. The summation (Sjj^fc) is over the SNR affected 
regions. Note that the change of mass will also affect the conservations of other values such as 

momentum and energy so this step must be done first. In general, fdens is very close to one, which 

is guaranteed by the dynamically generated SNR profile. 

(b) Conservation of momentum. There is no net momentum in each direction for the SNR 
seeds. However, the region to embed the SNR may contain non-zero momentum, especially when 
an SNR actually explodes in an outflowing wind region. Taking the momentum in the X direction 
as an example: the initial net momentum in X direction is 

MVXo = ^ Po (^> k)vx,o j, k)AVij^k , (A4) 
In order to conserver the momentum, we add a bulk motion 

Vxi = (A5) 
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to the embedded SNR velocity profile 

Vx,i{i,j,k) = Vx,sn{i,j,k) + Vx,i . (A6) 

MVXi = ^ P2 (i, J, k)vx,i (i, j, k)AVij,k , (A7) 

Here Vx,sn is the velocity profile of the SNR template. Let = MVXo /MVXi , thus 

Vx,2{i,j,k) = fvx X Vx,iii,j,k) . (A8) 

The same procedure is also applied to Y and Z direction. 

(c) Conservation of energy. The kinetic energy of the region is 

^-^0 = \ ^Po(^'^^^)['"x,oihj,k) +vl^o{i,j,k) +vlo{i,j,k)jAVij^k (A9) 
After adjusting the mass and momentum, the updated kinetic energy is: 

EK2 = ^Y.P^(^^^^'^H'"h(^^^^'^^ +vl2{iJ,kUAVij^k ■ (AlO) 

The energy conservation is achieved by adjusting the thermal energy accordingly: 

V2[%,3, k) = fsiPiiiJ, k), (All) 

where 

= ^^o + EK, + E^^-EK, ^ ^ ^ ^^°M^Ay.,,„ EI, = ^ ^^iM^AF,,,. (A12) 

Note the original random motion of mass inside the affected region is implicitly converted to thermal 
energy. The clumpy environment is also smoothed out by this scheme, and the embedded SNR is 
spherical symmetry related to the SNR center. 



